Critical viscoelastic response in jammed solids 
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We determine the linear viscoelastic response of jammed packings of athermal repulsive viscous 
spheres, a model for emulsions, wet foams, and soft colloidal suspensions. We numerically mea- 
sure the complex shear modulus, a fundamental characterization of the response, and demonstrate 
that low frequency response displays dynamic critical scaling near unjamming. Viscoelastic shear 
response is governed by the relaxational eigenmodes of a packing. We use scaling arguments to ex- 
plain the distribution of eigenrates, which develops a divergence at unjamming. We then derive the 
critical exponents characterizing response, including a vanishing shear modulus, diverging viscosity, 
and critical shear thinning regime. Finally, we demonstrate that macroscopic rheology is sensitive 
to details of the local viscous force law. By varying the ratio of normal and tangential damping 
coefficients, we identify and explain a qualitative difference between systems with strong and weak 
damping of sliding motion. When sliding is weakly damped there is no diverging time scale, no 
diverging viscosity, and no critical shear thinning regime. 
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Sufficiently dense packings of droplets in an emulsion, 
bubbles in a foam, or soft particles in a colloidal sus- 
pension form a jammed, mechanically rigid state [IH1]. 
Close to the (un)jamming point, which occurs at a crit- 
ical volume fraction <f> c 0.84 and 0.64 in d = 2 and 
3 dimensions 5-8,, contact deformations are small and 
particles are nearly spherical. Hence near unjamming 
these amorphous, viscoelastic solids can be modeled as 
athermal packings of soft spheres interacting via elastic 
and viscous contact forces |6], as in Fig. [I]}. 

One noted hallmark of the unjamming transition is the 
scaling of the quasistatic shear modulus Go, which van- 
ishes on approach to <j) c as Go/fc ~ (<j> — </> c ) 1 ^ 2 , when 
measured in units of the microscopic contact stiffness k 
[7] [59] . This anomalous scaling reflects the inherent soft- 
ness of materials close to unjamming and presages their 
loss of rigidity. Though unquestionably present in nu- 
merics |7, 9 13 , the scaling of Go has yet to be observed 
experimentally; in fact there have been few experimental 
results to date that unambiguously confirm scaling rela- 
tions emerging from the jamming paradigm [3J. This is at 
least partly due to the idealized nature of the frictionless 
soft sphere model; though granular matter is a common 
testbed for jamming [T4"HTrj] . friction qualitatively alters 
the way materials jam |17H19j . Foams and emulsions are 
closer realizations of frictionless soft spheres but present 
their own experimental challenges; quasistatic measure- 
ments, for example, are sensitive to coarsening and other 
sorts of aging on long time scales (2U] . 

The best opportunity to confront theoretical and ex- 
perimental studies of jammed viscoelastic solids may re- 
side in their finite rate response |12[ !21H32j , which is ex- 
perimentally accessible and physically significant in its 
own right [53Tl4"2"] . Linear viscoelastic response [T2J [S3] 
introduces a finite rate while keeping the strain ampli- 
tude infinitesimal; it bridges the gap between linear qua- 
sistatic response and response at finite strain and rate, 




FIG. 1: (a) Packing of soft repulsive spheres close to unjam- 
ming. Particles with no force-bearing contacts ("rattlers") 
have been removed. The mean coordination is z = 4.03, close 
to the isostatic bound z c — 4. (b) Relative motion of two 
disks in the center of mass frame of the upper disk. Parti- 
cles exert an elastic force proportional to their overlap Sij. 
Viscous forces oppose the relative normal velocity Aiiy and 
relative tangential velocity Aiijj — (pi8i+ pjOj) at the contact. 

including steady flows [Q l2lH24l 126] I29H3T] and shocks 
[32] . Here we measure numerically and calculate ana- 
lytically the complex shear modulus of jammed viscous 
sphere packings, 

G* (w) - G'H + iG"(uj) = 44 ■ (1) 

7(cj) 

G*(w) describes the linear response of a material sub- 
jected to oscillatory shear at frequency u), where ct{uj) is 
the stress amplitude and j(uj) is the strain amplitude. 
Its real and imaginary parts, G' and G", are respectively 
known as the storage and loss modulus, while the zero fre- 
quency limits of G' and G" /ui give the quasistatic shear 
modulus Go and dynamic viscosity i]q. 

This work builds on and substantially broadens results 
reported in Ref. [12]. We demonstrate that viscoelastic 
response near unjamming is fundamentally related to the 
relaxational modes of overdamped packings, and that the 
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FIG. 2: (a) Linear viscoelastic response is governed by two 
control parameters, the driving frequency cj and the excess 
packing fraction A<f> or, interchangeably, the excess coordina- 
tion Az or pressure p. There are three qualitatively distinct 
regimes of response distinguished by the degree of non-affinity 
of the particles' displacements (bl-b3; arrows indicate the di- 
rection and magnitude of a particle's motion under the pure 
shear depicted in b2), and the scaling of the complex shear 
modulus (c), shown here for the packing in Fig. [ij,. 



complex shear modulus is ultimately determined by the 
density of states D(s) of relaxational eigenrates s. We 
show that displacements at low frequency are consistent 
with quasistatic response when ujt* <C 1, with a crossover 
time scale 

* TO 



(2) 



that diverges at unjamming when sliding contacts 
strongly dissipate energy. The microscopic time scale 
tq is determined by coefficients in the viscous and elas- 
tic force laws. Because r* diverges at </> c , there is no 
quasistatic response at the unjamming transition. 

For times shorter than r* the dynamics enter a critical 
rate-dependent regime in which 

G*-fc(zwr ) A (3) 

with A = 1/2. Equivalent l/t A decay of the shear re- 
laxation modulus G r (t) has previously been observed in 
simulations of athermal suspensions close to jamming 
|25j . We refer to this as the critical or shear thinning 
regime, because doubling the driving rate increases the 
shear stress by less than a factor of two [50] . At high fre- 
quencies loto ^ 1, particle displacements smoothly follow 
the deformation gradient and the storage and loss moduli 
reflect the microscopic stiffness and damping coefficients. 
These quasistatic, critical, and affine regimes of response 
are summarized in Fig. [2j 

The anomalous quasistatic response near unjamming 
has its origins in the strongly non-affine nature of de- 
formations in jammed packings [8 . Non-affine motion 



occurs when individual particles do not smoothly follow 
the macroscopic deformation gradient; see Fig. [2}). Non- 
affinity manifests microscopically as an anomalous abun- 
dance of relative tangential motion at contacts, i.e. slid- 
ing, over relative normal motion H3] . The preference 
for sliding is closely linked to the character of floppy 
modes, zero energy deformations with no relative normal 
motion [5J H4Tl46| . While true floppy modes only exist 
in the unjammed phase, low energy excitations in the 
jammed state, which dominate quasistatic deformation, 
resemble floppy modes insofar as they have small normal 
motions and large sliding motions |47j . 

We will show that in viscoelastic response, non-affinity 
develops in time; the slower the deformation, the stronger 
the dominance of sliding motion. The three regimes of 
response are characterized by different degrees of non- 
affinity, as illustrated in Fig. [2^i,b. For a given packing 
fraction, non- affinity is maximal in the quasistatic limit, 
while in the affine regime sliding and normal motions 
are comparable. These two limits are connected by the 
critical regime, where non-affine motion is present but 
not fully developed. 

Because sliding plays a dominant role in low frequency 
response, the loss modulus is sensitive to the microscopic 
coupling between dissipation and sliding. We will show 
that a qualitatively different rheology results when slid- 
ing is weakly damped. Hence details of the microscopic 
interactions, in this case the local viscous force law, can 
dramatically alter the global rheology - in sharp contrast 
to the usual scenario for critical scaling in equilibrium 
systems [48] . 



I. LINEAR VISCOELASTIC RESPONSE OF 
SOFT SPHERES 

We begin by developing the equations of motion de- 
scribing linear viscoelastic response. Microscopic elas- 
tic and viscous interactions are linearized about a 
numerically-generated reference configuration, in this 
case a static soft sphere packing such as the one in Fig.[lji. 
To make numerical measurements, we average over en- 
sembles of packings prepared at different confining pres- 
sures p. Average properties such as the complex shear 
modulus G*(ui) and the relaxational density of states 
D(s) have characteristic features that depend only on 
the distance to jamming, e.g. the pressure p or the excess 
volume fraction A(f> := <j) — <f> c . We will determine these 
features using scaling arguments. While our numerical 
methods require computer-generated packings as input, 
the predicted scaling relations depend only on global pa- 
rameters such as the volume fraction or pressure. 



A. Equations of motion 

When a packing in d — 2 dimensions is deformed, the 
particles undergo displacements {ui}i=i...N from their 
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equilibrium positions. We consider viscous forces ca- 
pable of applying torques, which couple displacements 
and rotations {Oi}i=i...N ■ For macroscopic deformations, 
the packing also undergoes some pure shear strain 7 
(Fig. [2}d). The deformation is thus characterized by 
3N + 1 degrees of freedom: the vector components of the 
displacement, the particle rotations, and the strain 7. Al- 
though we present numerical results for two dimensions, 
all known critical exponents characterizing jamming are 



The dissipation function is [T2] 
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independent of dimension [8]. We will show in Sec 
that this is true for viscoelastic response, as well. 

We study soft particles interacting via repulsive elastic 
contact forces; this is the canonical model system for the 
study of (un)jamming 018]. The elastic force between 
disks of radii Pi and pj is linearly proportional to their 
overlap Sij = pi + pj — Ar^ , where Ary is the center-to- 
center distance. Overlapping disks also experience vis- 
cous forces that oppose relative normal and tangential 
motions at the contact; see Fig. [IJ>. We are interested 
in slowly driven systems, for which inertial effects can be 
neglected. Setting the particle masses to zero yields over- 
damped equations of motion. Overdamped soft spheres 
are often referred to as the "bubble model" after Durian, 
who introduced them as a model for wet foams [B]. 

The Lagrangian equations of motion for the bubble 
model are 



x ^ fdU &R _ 



(4) 



Here the {<9;}i=i...3A r +i label all the degrees of freedom 
of the system: the vector components of the particle dis- 
placements, the particle rotations, and the shear strain. 
U is the elastic potential energy and R is the Rayleigh 
dissipation function, both defined below. is the gener- 
alized force conjugate to the degree of freedom Qi. The 
generalized forces conjugate to the the displacements and 
rotations represent body forces and couples, respectively; 
they are all zero. The generalized force conjugate to the 
shear strain is the force moment aQ, where a is the shear 
stress and is the volume of the unit cell. 

We linearize the equations of motion Q about a static 
packing, while distinguishing the contributions of relative 
normal and tangential motions [44] . To quadratic order 
in 7 the elastic potential energy U is [61] 
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Here Au^ — (uj — vii 

disks along their contact normal hij , and Au 



fiij is the relative motion of two 



Au.,, 



Att|.- is the relative tangential motion; see Fig. lp. We 
take the contact stiffness kij = k to be the same for every 
contact. Quantities evaluated in the reference state are 
labeled with a superscript 0, and is the force carried 
by contact (ij) in the reference state. The contribution 
proportional to these forces is known as the pre-stress 
term. 
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One can verify that derivatives of R with respect to the 
particle velocities yield viscous forces damping normal 
and sliding motion at the contact with damping coeffi- 
cients 6" and respectively. Note that two overlapping 
but otherwise isolated particles separate with a relax- 
ation rate k/b" = 1/tq. In the following we shall express 
the damping coefficients in terms of the microscopic time 
scale To: fo" = kro and b 1 - = /3kro- The dimensionless 
quantity /3 is the ratio of tangential to normal damping. 

Collecting all degrees of freedom in a vector \Q(t)), 
the linearized equations of motion can be expressed as a 
matrix equation, 



t\Q(t))+B\Q(t))=a(t)\a), 



(7) 



where \a) = |{0},i7). The stiffness matrix K. and damp- 
ing matrix B capture the elastic and viscous interactions 
of the particles, 



d 2 U 
dQ m dQ ri 



and B mn = 



d 2 R 
dQ m dQ n 



(8) 



B. A note on units and scaling 

The stiffness matrix K. has units of the microscopic 
stiffness coefficient k, and the damping matrix B has 
units of the damping coefficient kro. In many physical 
systems, k and/or Zero depend on the amount of overlap 
5 between contacting particles [301 EHj • The typical over- 
lap itself vanishes at unjamming: it scales as S ~ A(j) 
[8]. Thus macroscopic observables can inherit a "triv- 
ial" dependence on the distance to unjamming directly 
from the microscopic coefficients. In the Hertz contact 
problem, for example, the pair potential V oc S a with 
a = 5/2, so the the contact stiffness k = V" oc 5 1 / 2 [49] . 
As noted above, G$/k ~ A0 1 / 2 near unjamming. Hence 
Go ~ Aifi in a system of particles interacting via Hertzian 
potentials, which is different from the square root scaling 
observed for harmonic interactions (a = 2). In our nu- 
merics, k and t are always set to unity; nevertheless, for 
clarity and generality, we explicitly indicate the k and To 
dependence of scaling relations. We omit dependence on 
the microscopic length scale p, the typical particle size, 
because it does not scale. 

As we are interested in scaling near unjamming, we 
must choose an appropriate measure of distance to the 
transition. We have already introduced A<p and p, which 
both vanish at unjamming; below we will also encounter 
the contact number z = z c + Az, which characterizes the 
geometry of a packing and approaches the isostatic value 
z c = 2d at unjamming [T3J [44] . Each measure has its 
own context-dependent advantages. We rescale our nu- 
merical measurements with Az because, as shall become 
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apparent below, it is the most fundamental of the three. 
However, as contact numbers are difficult to determine 
experimentally, it can also be convenient to couch scaling 
relations in terms of the dimensionless quantities A(f> or 
p/k. All three quantities are related via well-known scal- 
ing relations, Az 2 ~ Acf> ~ p/k [8], so scaling in terms of 
Az can always be rewritten as scaling in A<p or p/k. 



II. DYNAMIC CRITICAL SCALING 

We have shown that the viscoelastic shear response of 
a soft sphere packing obeys Eq. ^ . We will now use this 
equation of motion to determine the response to oscilla- 
tory driving at finite frequency. For a small oscillatory 
shear stress a = a(uj)e tult with frequency to, the particles 
oscillate with the same frequency, \Q(t)) = e lwt \Q(uj)). 
In general the response is out of phase with the stress. 
The equation of motion for oscillatory rheology in over- 
damped soft spheres is then 



JC\Q(u)) + iuB\Q(u)) = *(u)\&) 



(9) 



Eq. @ can be solved numerically for \Q(u)). Of par- 
ticular interest is the complex shear modulus, G* (w) = 
a(ui)/^/(u}), which provides a fundamental characteriza- 
tion of the shear response of the packing. The frequency 
dependent shear strain 7(0;) can simply be read off from 
|Q(w)), namely 7 (w) - Q-^alQiu)). 



A. Numerical results 

We evaluate the complex shear modulus numerically 
in an ensemble of packings generated with the molecular 
dynamics protocol of Somfai et al. [50]. Each packing 
contains N = 1024 weakly polydisperse particles prior to 
the removal of "rattlers," i.e. non-force bearing particles. 
For each excess contact number Az the response is aver- 
aged over ss 50 packings. Here and below, the pre-stress 
term of the elastic energy U has been neglected when 
evaluating the stiffness matrix IC; this approximation is 
justfied in Section |III| We find similar results when the 
pre-stress term is included, albeit with greater scatter. 

The complex shear modulus for varying excess coordi- 
nation Az is plotted in Fig.[3^i. To interpret the figure, it 
is helpful to recall that the simplest possible viscoelastic 
solid - the Kelvin- Voigt solid - has complex shear mod- 
ulus G^ v = Gkv + iVkvu, i.e. its storage modulus is flat 
and its loss modulus is linear in frequency. For both low 
and high frequencies, the packings of Fig. [3^, indeed have 
flat storage moduli. At high frequencies the storage mod- 
ulus is consistent with the microscopic spring constant, 
G'/k ~ 0(1)- The low frequency plateau in G' , however, 
does not match the high frequency plateau; the storage 
modulus is softer at low frequencies than at high, and 
therefore softer than one would naively anticipate based 



on the contact stiffness. This mismatch grows as the un- 
jamming transition is approached: the height of the low 
frequency plateau diminishes. 

The low and high frequency plateaus in the storage 
modulus are connected by a range of apparent power law 
scaling over a frequency interval of finite width. The 
power law is sublinear, and therefore the material is shear 
thinning. The upper bound of the shear thinning regime 
occurs for ujt ~ 0(1), independent of the distance to 
jamming. In contrast, the lower bound decreases - and 
the shear thinning regime grows wider - as unjamming 
is approached. 

Similar trends are observed in the loss modulus. Like 
the Kelvin- Voigt material, G" is linear at low and high 
frequencies. However, it also displays an anomalous 
power law regime for intermediate frequencies. The ap- 
parent exponent is less than one and similar to that ob- 
served in the storage modulus. The upper and lower 
bounds of the shear thinning regime coincide with those 
observed in the storage modulus. 



B. Scaling 

We noted above that the quasistatic (u> — > 0) shear 
modulus Go of soft sphere packings scales as 



Go 

T 



Aft 



Az 



(10) 



with \i R3 1.0 P3Q2]. Thus G vanishes continuously at 
the unjamming transition. In many soft matter systems 
with a rigidity transition, finite frequency response near 
the transition displays dynamic critical scaling [48l [51] . 
We now show that this is also true of viscoelastic solids 
near unjamming. First we demonstrate scaling numeri- 



cally; then, in Section III we explain its origins. 
We begin by making the scaling ansatz 



i G .(„,A=)=A,.. S .(0) 



(11) 



where the scaling function G*(x) = G'(x) +iG"(x) obeys 



and 



G'(x) 



Q"{x) 



const i<Cl 



x 
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1. 



X I«l 
x A X > 1 



(12) 



(13) 



The x <C 1 form of the dimensionless scaling function 
G* is chosen to reflect the frequency dependence of the 
Kelvin- Voigt material. The same exponent A appears in 
both G' and G" for x ^> 1; this follows from the Kramers- 
Kronig relations [62]. To guarantee that G*(uj > 0) re- 
mains finite when Az = 0, the exponents /i, A, and A 
must obey 



/i = AA . 



(14) 
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FIG. 3: (a) Storage modulus G' (open symbols) and loss modulus G" (filled symbols) for a range of excess coordination 
Az (legend) close to unjamming. Normal and tangential damping coefficients are equal (/3 — 1). The dashed line has slope 
1. (b) Critical collapse of the low frequency data in (a). The short dashed line has slope 1; the long dashed line has slope 
A = fi/\ = 0.54. 



Dynamical critical scaling relates two qualitatively dif- 
ferent regimes of response. The rheology of Fig. [3^,, 
however, clearly displays three regimes - these are the 
quasistatic, critical, and affinc regimes of Fig. [2) In the 
following we seek to relate the quasistatic and critical 
regimes; we therefore exclude frequencies ujt > 0.3. To 
test the critical scaling ansatz, we plot (G* /k)/ Az^ 1 ver- 
sus (ijjto)/Az x and vary fi and A to find the best data 
collapse. We find /i = 1.05(5) and a dynamical exponent 
A = 1.95(5). The resulting collapse is excellent, as shown 
in Fig.[3j3. The shear thinning exponent A = fi/X = 0.54 
is also in good agreement with the data. In Section [TT] we 
will predict /j, = 1, A = 2, and A = 1/2, but for now they 
are phenomenological. 

The data collapse in Fig. ^jp indicates that low fre- 
quency dynamics near jamming display dynamic critical 
scaling. This has several immediate implications. First, 
the scaling function has a characteristic time scale 



The quasistatic limit corresponds to frequencies u <C 
1/t*. In this limit we can read off the quasistatic shear 
modulus Go and dynamic viscosity rjo- The shear mod- 
ulus is Go = lim w _j.o G'(w) ~ kAz^, which recovers the 
quasistatic scaling of Eq. ( p^Oj ) as expected. The dynamic 
viscosity is determined by the imaginary part of the scal- 
ing function, 770 = lfin^^o G"(u))/uj. It scales as 

% ~ G r* ~ , (16) 

which is diverging. Note that the viscosity does not di- 
verge with the same exponent as the time scale r* due 
to the vanishing shear modulus. 

The storage and loss moduli are dominated by the qua- 
sistatic moduli for frequencies up to 1 jr* . On time scales 
shorter than r*, however, the dynamics enter a qualita- 
tively different, critical regime. Instead, both the storage 



and loss moduli display an anomalous power law scaling 
with frequency, expressed succinctly as 

G*-fc(iwT ) A . (17) 

Because < A < 1, this response is a form of shear 
thinning. A is a new critical exponent near unjamming. 

Our low frequency scaling relations can be compared to 
known results calculated in the strictly quasistatc limit, 
which corresponds to setting lo to zero in Eq. {9). If 
the limit is continuous, then quasistatic particle trajec- 
tories should be a good approximation to response at 
sufficiently low frequencies, where damping is weak. We 
noted above that the exponent fj, characterizing the low 
frequency shear modulus is in good agreement with the 
established quasistatic value of 1.0, which suggests the 
quasistatic limit is indeed continuous. Under quasistatic 
deformation the typical magnitude of sliding motion is 
known to scale as Am 1 - ~ 7/Az 1 / 2 [5] S3]. This can 
be related to the viscosity by balancing the dissipated 
power (l/2)?7o(w7) 2 with the microscopic dissipation rate 
/ visc -A-it - fcrofwAu 1 ) 2 , giving 770 ~ kr /Az. Note that 
the quasistatic regime breaks down when t]qui/Go <C 1 
is violated, i.e. when the quasistatic scaling relations no 
longer permit damping to be weak compared to storage. 

III. RELAXATIONS AND RESPONSE 

In Section [TT] we demonstrated numerically that shear 
response near jamming displays dynamic critical scaling. 
Simply demonstrating scaling does not explain its origin, 
so we now seek to identify the physical mechanism under- 
lying packings' critical response. We begin by showing 
that viscoelastic response is governed by the relaxational 
cigenmodes of an overdamped packing. We then numer- 
ically determine the relaxational density of states D(s) 
and introduce a scaling argument to explain the charac- 
teristic features of D(s). Finally, we demonstrate that 



6 



the exponents fi, A, and A characterizing oscillatory rhe- 
ology can be derived from the density of states. 



Sn-l 

1 



s n+l 



A. Complex compliance and shear modulus 

Anticipating that the eigenmodes of an overdamped 
packing are evanescent, it is natural to consider the 
Laplace transform of the equation of motion. In the ab- 
sence of driving, a = 0, this is 



(lC + s n B)\s n ) =0. 



(18) 



Eq. ( 18 1 is a generalized eigenvalue equation - "general- 



ized" because B is not the identity matrix. The equa- 
tion has eigenvalues {s n } and orthonormal eigenvectors 

Eq. ( |18[ ) is invariant under rigid body translations and 
individual particle rotations; thus there are N + 2 trivial 
zero modes of the system in two dimensions. A jammed 
system does not possess non-trivial or "floppy" zero 
modes. The remaining eigenvectors are 23-orthogonal, 
meaning {s m \B\s n ) — for m ^ n, and the remaining 
eigenvalues are negative definite. Physically, the latter 
property means that the modes are evanescent: a pack- 
ing deformed along a mode |s„) relaxes exponentially to 
the reference configuration with relaxation rate \s n \. For 
convenience, from this point forward we adopt a sign con- 
vention in which the symbol s n refers to the absolute 
value of the eigenrate of the n th eigenmode. 

We now return to Eq. ([£]) with driving. Expressing the 
response as a superposition of relaxational eigenmodes, 
one may invert the equation of motion to solve for the 
complex compliance J*(oj) = j(u>)/a(u>), 



n ' 



T) n LU 2 + S 2 n 



r) n ui 2 + si 



(19) 



We have introduced r}„ = k(s n \B\s n ) /\(s n \a)\ 2 ; the mo- 
tivation for this nomenclature will become apparent be- 
low. Note that the projection (s n \a) is a measure of how 
strongly mode n couples to shear. The trivial transla- 
tion^ and rotational modes do not couple to shear. 
Eq. ( 19 1 completely specifies the microscopic response 



of the system in terms of properties of the individual 
modes, including the complex shear modulus G*(w) = 
l/J*(ui). It is straightforward to show that response of 
this form is identical to that of a series circuit of viscoelas- 
tic elements, one for each mode, as depicted in Fig. [4] 
Such a circuit has a complex shear modulus 



G* = 



^G* 



(20) 



Each element has a modulus G* = G n +ir) n uj, where r/ n is 
defined above and G n = i)„s n . Hence each Kelvin- Voigt 



FIG. 4: Series circuit of Kelvin- Voigt elements. Each element 
has a characteristic relaxation rate associated with one of a 
packing's relaxational modes. 



element has a characteristic relaxation rate G n /rj n = s n 
equal to the eigenrate of the corresponding mode. We 
stress that the series circuit of Fig. [4] is not an ad hoc 
model, but a mapping from an exact solution to the equa- 
tions of motion. 

We have now shown that a packing's complex shear 
modulus is determined by its relaxational modes, and 
that two quantities are required to quantify each mode: 
its relaxation rate s n and the coupling of the n th mode to 
shear via the mode's "viscosity" r) n . To make progress, 
we now introduce an approximation: we assume that r\ n 
is random, i.e. independent of n and implicitly of s n as 
well, and that it has a typical value r\ n ~ kr^. We have 
verified numerically that this is reasonable. 

Under the above assumption, the fact that response 
can be mapped to the series circuit of Fig. [4] establishes 
an important intuition. The modes with the lowest eigen- 
rates, i.e. those that relax slowly, have the smallest effec- 
tive storage modulus G n ~ kros n . As a series circuit's 
low frequency deformation is carried by its softest ele- 
ments, one infers that low frequency shear will be carried 
predominantly by the slow modes. 



Taking the thermodynamic limit of Eq. ( 19 1 gives 



kf(oj) 
fcj"(w) 



1 



TO Jo « 



T JO 



^D(s)ds 

'—D{s)ds. (21) 



The scaling of J*, and hence G* as well, is now com- 
pletely governed by the relaxational density of states 
D(s), which we proceed to determine. 



B. Relaxational modes 

What are the characteristic features of D(s)7 And how 
are they reflected in the response? To answer these ques- 
tions, we investigate D(s) both numerically and theoret- 
ically. We will show that slow modes are both abundant 
and strongly non-affine. 

Fig. [5^i depicts the relaxational density of states for 
several different distances to jamming, parameterized by 
Az. Several features stand out. First, the fastest relax- 
ations are on the order of the bare rate 1/tq at which 
two overlapping particles separate. Second, low-s relax- 
ations carry a large statistical weight, and that weight 
increases as jamming is approached: the closer it is to 
jamming, the more ways a packing has to relax slowly. 
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FIG. 5: (a) Ensemble averaged relaxational density of states D(s) near the unjamming transition, where the excess coordination 
Az (legend) vanishes. Normal and tangential damping coefficients are equal (/3 = 1). The dashed line has slope 1/2. (b) Rescaled 
data from (a). The exponents A and A have been chosen so as to collapse the crossover at low s and generate a plateau above 
it. 



Third, and more quantitatively, D(s) is approaching a 
power law divergence in the limit Az — > 0, 



TO 



1 

ST 



A' 



(22) 



for some exponent A' > 0. At finite Az the divergence 



A: 



for some 



is cut off at a characteristic rate s* 
A' > 0. 

To characterize the relaxational density of states, we 
plot s A D(s) versus s/Az x . We vary A' to find the value 
for which D(s) approaches a plateau and vary A' to find 
the value that collapses the low-s crossovers. We find a 
flat plateau and collapsing crossovers for A' = 0.47(5) 
and A' = 1.9(1). Results are shown in Fig. ^jp. Note 
that A' and A' are in good agreement with the critical 
exponents A and A, respectively, that govern response. 
Below we will show that A' = A and A' = A. 

To better understand the nature of the abundant low-s 
modes in the relaxational density of states, we relate their 
eigenrate to the non-affine character of their motion. The 
prevalence of sliding over normal relative motion can be 
quantified by the ratio 



(Am-L) 2 



(Aull 



(23) 



which is evaluated for the relative displacements of mode 
n. Bars indicate averages over all contacts. We now show 
that r„ becomes larger as s ra is decreased, indicating the 
growing importance of sliding motion. 
Each eigenrate of the material satisfies 



(3n\£\Sn) = U-U C 

(s n \B\s n ) 



R 



(24) 



where U and R are the elastic potential energy and 
Rayleigh dissipation function of Eqs. ([5| and ([6]), respec- 
tively. By considering typical values of the relative dis- 



scripts, the relation is 



placements, Eq. (24 1 relates s n to r n . Dropping sub- 

(25) 



ST 



T- 2 + c(p/k) 



dp 



where the prefactors c and d are constants expected to be 
of order unity. To reach Eq. (251, we have used the fact 



that the typical force in the reference state is proportional 
to the pressure and have assumed that Au 1 - dominates 
sliding at the contact. Specializing to the case j3 = 1, 
where normal and tangential motion are damped equally, 
and solving for T in the limit of small s gives 



1 

f2 



const 



stq 
(p/k) 



(26) 



Thus the slower the mode, the stronger the dominance of 
sliding over normal motion. Eq. ( 26 1 is verified in Fig 



which plots r for six packings (no ensemble averagin. 
at different confining pressures. The pre-stress introduces 
an upper bound on the non-affinity; T approaches a con- 
stant for relaxation rates s < (p/k)/ro, which coincides 
with s*. At the crossover rate s* sliding exceeds normal 
motion by a factor T* ~ 1/Az. 



C. Predicting the relaxational density of states 

The exponents A' and A' characterize the form of the 
density of states. We now seek to determine them ana- 
lytically. Several years ago, Wyart, Nagel, and Witten 
(WNW) used a variational method [47] to explain the 
anomalous plateau in the vibrational density of states 
D(Cj) - the probability density of vibrational eigenfrc- 
quencies Q of packings with inertia and without damp- 
ing. For the case where all particles have unit mass, the 
square of each eigenfrequency, w 2 , is an eigenvalue of the 
stiffness matrix K,. The undamped vibrational modes are 
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FIG. 6: The ratio Y = Au 1 /Au" as a function of relaxational 
eigenrate s collapses when rescaled with confining pressure p 
according to Eq. (26 1. The dashed line has slope one. 



therefore a natural basis in which to express quasistatic 
deformations. 

The work of WNW has become a landmark of the jam- 
ming literature. Here we generalize the method to over- 
damped dynamics and show that it correctly predicts the 
values of A' and A'. Without going into technical details, 
we sketch the key ideas of the method and show how they 
can be extended to overdamped dynamics. 

The central idea of the WNW method is as follows. 
One takes a packing of particles with contact number 
z — z c + Az and lays a mesh over it. The packing has 
linear size L and the mesh has a characteristic length 
q" 1 . Everywhere a contact crosses a face of the mesh, 
that contact is "cut" or ignored. If enough contacts are 
cut, the packing loses rigidity and zero energy collective 
motions, i.e. floppy modes, appear. For this to happen 
the number of cut contacts, which is on the order of the 
total surface are of the mesh O(qil), must exceed the 
number of contacts the original packing had in excess of 
the isostatic value, which is O(Azfl). Therefore floppy 
modes appear whenever q > q* ~ Az. The diverging 
length scale 1 jq* is known as the isostatic length (45J |47j . 
When q > q* the cutting process creates on the order of 
(q — q*)Q floppy modes. 

The floppy modes can be exploited to construct trial 
modes {|Qwnw)} for use in a variational argument. Each 
trial mode is a floppy mode modulated by a sinusoidal 
envelope that vanishes on the faces of the mesh. This "re- 
pairs" the large relative motions between particles par- 
ticipating in cut contacts. 

Having determined how to construct trial modes, the 
vibrational density of states can be determined in two 
steps. Noting that trial modes are parameterized by 
their "wavenumber" q, one first determines the proba- 
bility density D{q). The integral of D{q) is related to 
the number density of trial modes, J 9 , dq' D(q') ~ q — q*- 
Differentiating with respect to q gives 

D{q) ~ const for q > q* . (27) 

The second step is to relate the trial mode wavenumber q 



to the frequency Cj via a dispersion relation u = Cj(q); one 
then has D(w) = D(q)\dq/dQ\ ~ |dq/dw|. For undamped 
modes the dispersion relation turns out to be linear in q, 
so D(£j) develops a plateau. 

For our purposes, the key observation is that the WNW 
argument is exclusively geometric in nature, up to the 
point where the dispersion relation u)(q) is invoked. The 
undamped nature of the dynamics enters only in this fi- 
nal step. Recall that the undamped vibrational modes 
are the natural basis with which to describe quasistatic 
deformations. We anticipate that trial modes that work 
well for quasistatic deformations will also work well for 
sufficiently slow deformations. Therefore we can use the 
same WNW trial modes, parameterized by q, and ap- 
ply them to overdamped dynamics to estimate the relax- 
ational density of states, 

D(s)=D(q)\dq/ds\~\dq/ds\ for s > s* . (28) 



Here we have introduced the crossover rate s* = s(q*), 
which we identify with 1 /r* . 

It remains only to determine the overdamped disper- 
sion relation s = s(q). To do this we exploit Eq. (251, 



the relation between relaxation rates and non-affinity 
Recall that trial modes are constructed by applying a 
sinusoidal modulation to a floppy mode created by the 
cutting procedure. Because of the modulation, the rel- 
ative normal and tangential motions of the trial mode 
differ from those in the floppy mode. As the sinusoid has 
wavenumber q, its contribution to the relative motions is 
of order quf m , where Uf m is the typical displacement in 
the floppy mode. This amounts to a small correction to 
the tangential motions, which are nonzero in the floppy 
mode. However, it sets the scale of the relative normal 
motions of the trial mode, because the relative normal 
motions of the floppy mode are zero, to leading order. 
Therefore r WNW = Au^r NW /Ait| VNW ~ 1/q. 
The dispersion relation for trial modes is thus 

1 + c(p/k) 



st 



c'/3 



(29) 



Recalling that q > q* for all trial modes, the disper- 
sion relation can be expanded, keeping only the dominant 
term: 



ST 



q 2 /f3 
const 



q « P 1/2 

/3 1 / 2 «g. 



(30) 



As the the dispersion relation is independent of p to lead- 
ing order, the pre-stress can be neglected. When normal 
and tangential motions are damped similarly, /3 is of or- 
der unity and the dispersion relation is quadratic. That 
overdamped dynamics should display a quadratic disper- 
sion relation is not surprising; consider, e.g., plane waves 
in an overdamped ball-and-spring chain. We stress, how- 
ever, that this quadratic form was not assumed, but de- 
rived from properties of the trial modes. Indeed, when 
the damping ratio (3 < (<7*) 2 i i-e. when sliding motion 
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is weakly damped, q is always larger than /3 1 / 2 , the dis- 
persion relation is no longer quadratic, and there is a 
crossover to qualitatively different response. We return 
to this observation in Section IV. 

For later convenience we write the q -C /3 1 / 2 dispersion 
relation as 



ST 



Q 



for 



> 



(31) 



Proceeding in the same way, one finds that the dy- 
namic viscosity scales as rjo ~ Gqt* ~ kro/Az x with 
exponent A' — /i = 1. Hence A = A', as anticipated, and 
we have also explained the numerical finding of Eq. (16 1. 

Shear thinning regime. — For frequencies s* < uj < 
1/t , quasistatic predictions fail and the response be- 
comes shear thinning. The real and imaginary parts of 
the complex compliance are 



Hence X' = 2 and 



S T 



(<f) 



Y 



P 



Az x ' 



Invoking Eq. ( 28 1 , the density of states scales as 

A' 



D(s) 
Pr 



1 

/3sT 



for 



> 



(32) 



(33) 



with A' = 1/A' = 1/2. The values of both A' and A' 
are in excellent agreement with the numerical results of 



Section HI B| for /? = 1. We will confirm the scaling with 
P in Section |IV| 



kf 



and 



k J" 



1 

lA 7 



,1-A' 



ds 





A' 


LOTo J 




1 


' r ds 


w 




(± 







I/to 



,1+A' 



uj ds 

c2+A' 



(35) 



(36) 



D. Predicting the complex shear modulus 

With the form of the relaxational density of states in 
hand, we can now return to the complex shear modulus. 
We restrict our attention to the case = 1 until indicated 
otherwise. We shall show that the divergence in D(s) 
and the crossover rate s* , through the exponents A' and 
A', suffice to predict the scaling of G*(ui) demonstrated 
above. Specifically, we will show that /i = A' A', A = A', 
and A = A'. 



Eq. (21 1 relates the complex compliance to an integral 
over modes. Using our results for the relaxational den- 
sity of states, the integrals of Eq. (21 ) can be made into 
bounds on the scaling of J* by replacing J Q (-)D(s)ds 

with t J*T( - )( ST o) _A ds. We shall assume that the 
bounds are saturated, which is verified by the good agree- 
ment between the resulting predictions and the numerical 
results presented above. 

Quasistatic response. — We begin with the zero fre- 
quency or quasistatic limit. The quasistatic shear 
modulus Go and dynamic viscosity r/o of a viscoelas- 
tic solid obey 1/G = lim^^p J' {to) and rjo/G^ — 
♦o[J»/w] 



k 

Go 



Eqs. (|2lJ), (|31j), and (|33j give G : 



1 



«1+A' 



(q*) 



A' A' 



(34) 



Using q* ~ Az, it immediately follows that the qua- 
sistatic shear modulus scales as Go ~ k Az' 1 with \x = 
A' A' = 1. This explains the numerical finding of Eq. ( [To| ) 
and is in good agreement with prior numerics [3 19] . Our 
result is compatible with previous calculations of the qua- 
sistatic shear modulus [TTJ [55] and is, to our knowledge, 
the first to directly relate Go to the isostatic length 1 jq* . 



It immediately follows that G* ~ k(iujTo) A with A = 
A' = 1/2, again in excellent agreement with numerics. 

In linear response the stress relaxation modulus G r (t), 
which gives the stress response after a unit step strain at 
time t = 0, can be calculated directly from the complex 
shear modulus. The scaling G* ~ («w) A implies a regime 
of power law relaxation G r (t) ~ t~ A . Therefore our re- 
sults explain the critical relaxation G r ~ t~ 5 seen in 
simulations of athermal suspensions near jamming [25 . 

High frequency response. — The density of states van- 
ishes above an upper bound on the order of 1/tq. For 
driving frequencies that exceed the bare rate, one finds 
J' ~ l/k(uiTo) 2 and J" ~ — 1/kuJTQ. Hence the high fre- 
quency moduli are simply given by the coefficients of the 
microscopic force law, 



G' ~ k and G" - kujT . 



(37) 



This is consistent with the observation that high fre- 
quency response is affinc. 

To summarize, when expressed in terms of A</>, the 
storage modulus scales as 



G' 



k A^/2 

k^To) 1 ' 2 
k 



LjJTo < A(j) 

A<t> <lot <l (38) 

1 < UTo , 



and the loss modulus scales as 



G" 



k(ujT Q )/A^/ 2 

fc(wT ) 



LUTo < Alp 

A(f> <lot <1 (39) 

1 < LOTo ■ 
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FIG. 7: (a) Relaxational density of states for Az = 0.063 
and varying damping ratio j3 (legend) . (b) Evolution of the 
relaxational density of states for Az = 0.063 and varying /3. 
Inset: Dynamic (to — > 0) viscosity 770 and affine (cj — > 00) 
viscosity r/oo as a function of /3. 



IV. WEAKLY DAMPED SLIDING 

In the preceding Sections the damping ratio j3 was 
taken to be unity, so relative normal and tangential mo- 
tions were damped with equal strength. The dispersion 
relation, Eq. ( 30 ) , then has a simple quadratic form. This 



is significant because the form of the relaxational den- 
sity of states and the complex shear modulus both follow 
from the dispersion relation. As the dispersion relation 
ceases to be quadratic when sliding is weakly damped, 
one anticipates changes in D(s) and G*(uj). We there- 
fore investigate the limit of weak tangential damping. 

The ft — > limit is important for an additional rea- 
son. While it may seem obvious that sliding motions are 
damped in materials such as foams and emulsions, many 
numerical studies take /3 to be zero, so that only relative 
normal motions dissipate energy [23l [24j [26l [53] . This 
is done for reasons of numerical convenience: with slid- 
ing undamped, neither elastic nor viscous forces impose 
torques, so the particles' rotational degrees of freedom 



can be ignored. We will show that this choice of force 
law qualitatively alters a system's viscoelastic response. 

The condition for weakly damped sliding is /3 Az 2 . 
When sliding motion is weakly damped, the dispersion 
relation stq ~ const is independent of q to leading order 



- see Eq. (30). This already suggests a dramatic depar- 
ture from the dynamic critical scaling described above, 
as sending q — > q* no longer produces a diverging time 
scale. 

To see how response evolves as a system passes from 
strong to weak tangential damping, in Fig. [7^, we plot 
the relaxational density of states for varying j3 averaged 
over 20 packings with excess coordination Az = 0.063. 
Several features are apparent. For (3 > 1 there is a gap 
between slow modes with /3sr < 0(1), which has the 
square root divergence discussed above, and a narrow 
band with stq ~ O(l). D(s) narrows with decreasing 
/3 < 1, consistent with the form of the dispersion rela- 
tion. This is because the crossover rate s* ~ Az 2 / (3tq 
increases, and in so doing "squeezes out" the anoma- 
lous slow modes: the interval over which D(s) displays 
l/s A growth becomes increasingly narrow. In the inset 
of Fig. [7}l we plot D/(3t versus /3stq for slow modes 
in the strong tangential damping regime (sro < 0.1 and 
(3 > Az 2 ). The data collapse is excellent, confirming the 
j3 dependence of Eq. (33 1. 

In Fig. (JtJd) we plot the complex shear modulus for 
varying (3, Clearly the form of G* is dramatically differ- 
ent when sliding is weakly damped. In particular, the 
critical shear thinning regime narrows as j3 — ¥ and 
vanishes when s*tq ~ O(l), which coincides with the 
weak tangential damping regime. Both the dynamic vis- 
cosity r/Q ~ j3 and the high frequency (affine) viscosity 
Voo ~ /3 1 / 2 depend on the damping ratio when sliding 
is strongly damped (Fig. [7]b inset) , but become indepen- 
dent of (3 in the weak tangential damping regime. To 
understand the crossover in 770, recall that when slid- 
ing is strongly damped, the dynamic viscosity 770 ~ 
(Au J -/7) 2 ~ (3/Az reflects the scaling of the typical 
relative tangential displacement in quasistatic response. 
When sliding is weakly damped, the dominant source of 
dissipation is relative normal motion and ryo ~ (Au" /7) 2 , 
independent of f3. Ellenbroek et al. find Am" /j ~ Az 1 / 2 
in the quasistatic limit [9]. 

To gain further insight into the weak tangential damp- 
ing regime, it is useful to directly consider the case f3 — 
while also neglecting the pre-stress. The pre-stress can 
be neglected because it does not appear in the leading 
order term in the dispersion relation. The key observa- 
tion is that under these conditions, the damping matrix 
B and the stiffness matrix K, are proportional, B = tqK,. 
The equation of motion becomes 

[l+UJT )fc\Q{u)) = a{u)\&). (40) 

By expanding \Q(w)) in th e un damped eigenmodes of 
the stiffness matrix IC, Eq. (401 can be solved for the 



complex shear modulus. The result is 
G*(oj) = G (1 + iojt ) 



(41) 
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The quasistatic modulus Go ~ k Az is unchanged, but 
the dynamic viscosity ryo = GqTq now vanishes, as antici- 
pated from the quasistatic scaling of Am" . The time scale 
Vo/Go, which equals r* when sliding is strongly damped, 
no longer diverges. Moreover, consistent with our nu- 
merical observations, there is no critical shear thinning 
regime. 

The above discussion establishes an important intu- 
ition. Quasistatic response near unjamming is strongly 
non-affine [43] , reflected in the critical divergence of 
sliding motion on approach to the unjamming point. 
When non-affine motion dissipates energy, its divergence 
is reflected in the viscosity; see also Refs. [12 02HSS]- 
Dynamic critical scaling, which describes response at fi- 
nite rates, requires a diverging time scale. We have 
shown that non-affine motion must be strongly damped 
to have a diverging time scale in viscoelastic linear re- 
sponse. When sliding is damped weakly or not at all, 
there is no diverging time scale and hence no dynamic 
critical response, reflected in the absence of a shear thin- 
ning regime. 

V. CONCLUSIONS 

Linear response in viscoelastic jammed solids displays 
dynamic critical scaling when sliding is strongly damped. 
We have shown that critical response is driven by slow 
relaxational cigenmodes. Close to unjamming the slow 
modes are anomalously abundant and have a spatially 
non-affine character reminiscent of floppy modes. By ex- 
tending the WNW method to overdamped dynamics, the 
broad distribution of relaxation rates can be explained. 
The complex shear modulus can then be derived from 
the relaxational density of states. 

Our results clearly establish that nonlinear viscous 
forces are not necessary for a material to be shear thin- 
ning and, concomitantly, that low frequency rheology 
near jamming does not trivially reflect the form of the 



microscopic viscous force law. We have also shown that 
viscoelastic response displays a sensitive dependence on 
the form of the viscous force law. The critical shear thin- 
ning regime only exists when sliding motion is strongly 
damped, and in the weak tangential damping regime the 
dynamic viscosity vanishes on approach to unjamming. 
We consider it likely that the qualitative distinction be- 
tween strong and weak tangential damping persists in 
nonlinear response, including steady flow. Therefore sim- 
ulations that damp only normal motion EH] EH ESI [53] 
may not be comparable with those that include damping 
of tangential motion [2UH31j . 

Our calculations can be compared with a number of ex- 
perimental results [33l[35l|36llMl[39lll2]- As noted above, 
the quasistatic regime has yet to be observed. Foams and 
emulsions, for example, are not thermodynamically sta- 
ble; their structure evolves due to coarsening and related 
effects not included in the bubble model. As a result, 
their low-frequency rheology depends on the age of the 
sample ED]- More encouragingly, there is experimental 
evidence of oj - 5 scaling in the complex shear modulus of 
emulsions ; 33] , liquid foams [35J [SHI I3S] , organic foams 
[55] . and soft suspensions consistent with our pre- 
diction for G* in the critical regime. This comparison 
must be made with caution, however, because the vis- 
cous force law is not always known and likely varies from 
material to material. While viscous forces are likely to be 
linear in some of these materials, the bubble-bubble drag 
force in foams is known to be nonlinear (3TJ [57] [58] . An 
important goal for future research will be to identify the 
full dependence of macroscopic rheology near jamming 
on the form of the microscopic force laws, including non- 
linear viscous forces. 
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The stiffness kij of the contact between particles i and j is 
the second derivative of the elastic pair potential V(rij), 
where rij is the distance between the particles' centers. 
Conventionally the term "shear thinning" is reserved for 
steady flow, but we are not aware of a suitable analog for 
oscillatory rheology. 

In general, Eq. |5| also has a term a°Q^y that is linear in 
the strain 7. However, we consider only isotropic refer- 
ence states, for which the shear stress <r° = 0. 
To see this quickly, consider the local Kramers-Kronig 
relation G" ~ (ji/2)u) d^G' , which is appropriate when 
the storage modulus is weakly dependent on frequency. 
More careful consideration of the full (integral, hence 
non-local) Kramers-Kronig relations leads to the same 
conclusion. 



